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^sj Abstract. The elliptic Monge-Ampcre equation is a fully nonlinear Partial 

Differential Equation that originated in geometric surface theory and has been 

applied in dynamic meteorology, elasticity, geometric optics, image processing 

^J and image registration. Solutions can be singular, in which case standard 

numerical approaches fail. Novel solution methods are required for stability 
r/*\ and convergence to the weak (viscosity) solution. 

In this article we build a wide stencil finite difference discretization for the 
I I Monge-Ampcre equation. The scheme is monotone, so the Barlcs-Souganidis 

^T theory allows us to prove that the solution of the scheme converges to the 

h^ unique viscosity solution of the equation. 

^H Solutions of the scheme are found using a damped Newton's method. We 

pH prove convergence of Newton's method and provide a systematic method to 

■^-^ determine a starting point for the Newton iteration. 

Computational results are presented in two and three dimensions, which 
demonstrates the speed and accuracy of the method on a number of exact 
I I solutions, which range in regularity from smooth to non-diffcrcntiable. 
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y *«0 1. Introduction 

^^ In this article we introduce a monotone discretization of the Monge-Ampere 

r — , equation, which is valid in arbitrary dimensions. A proof of convergence to the 

viscosity solution is provided, as well as a proof of convergence of Newton's method. 

Numerical results are also presented. 
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1.1. The setting for equation. The Monge-Ampcre equation is a fully nonlinear 
Partial Differential Equation (PDE). 



(MA) det(D 2 u(x)) = fix), for x in ft. 

03 

The Monge-Ampere operator, det(D 2 u), is the determinant of the Hessian of the 
function u. The equation is augmented by the convexity constraint 

(C) u is convex, 

which is necessary for the equation to be elliptic. The convexity constraint is made 
explicit for emphasis; it is necessary for uniqueness of solutions and it is essential 
for numerical stability. 
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While other boundary conditions appear naturally in applications, we consider 
the simplest boundary conditions: the Dirichlet problem in a convex bounded sub- 
set SlcK 1 ' with boundary dil, 

(D) u(x) = g(x), for x on dil. 

Under suitable assumptions on the domain and the functions f{x),g(x), recalled 



subsection 2.1 



there exists a unique classical (C 2 ) solution to (MA}, dC), (|d|) 



However, when these conditions fail, the solution can be singular. For singular 
solutions, the correct notion of weak solution must be used. In this case, novel 
discretizations and solution methods must be used to approximate the solution. 



1.2. Applications. The PDE (MA) is a geometric equation, which goes back to 
Monge and Ampere (see [16). The equation naturally arises in geometric problems: 
the Minkowski problem for prescribing the Gaussian curvature of a surface [51 H0l[3T)] 
and the related isometric embedding problem for manifolds 29J . Early applications 
identified in [39] include dynamic meteorology, elasticity, and geometric optics. 

The Monge- Ampere equation arises as the optimality condition for the problem 
of optimal mass transport with quadratic cost [TH1 HI US] . This application of the 
Monge- Ampere equation has been used in many areas: image registration [27J I2H1 
El3"] , mesh generation [T51 ITUl [5] , reflector design [53] , and astrophysics (estimating 
the shape of the early universe) [20] . 

The optimal transportation problem seeks a map g(x) that moves the measure 
m( x ) to ^2(2/) arL d minimizes the transportation cost functional 

x-g(x)\ dm. 



The optimal mapping is given by g = Vw, where u is convex and satisfies the 
Monge- Ampere equation 

dct(D 2 u(x)) = /ii(a;)//i 2 (Vii(i)). 

In this situation, the Dirichlet boundary condition |D| is replaced by the implicit 
boundary condition 

(1) g(-):"l->n2, 

where the sets fii s O2 are the support of the measures /ii,/i2, respectively. These 
boundary conditions are implicit and thus difficult to implement numerically. For 
many applications, both domains are squares, and a simplifying assumption that 
edges are mapped to edges allows Neumann boundary conditions to be used. In 
other applications, periodic boundary conditions are used. Building on the work 
in this article, the first author has recently devised a method to implement the 
boundary condition (fTl) [22] , 

In other problems, the Monge-Ampere operator appears in an inequality con- 
straint in a variational problem for optimal mappings where the cost is no longer 
the transportation cost. Here the operator has the effect of restricting the local 
area change on the set of admissible mappings; see [53] or [TT] . 

1.3. Related numerical works. Despite the number of applications, few pub- 
lications devoted to the numerical solution of the Monge-Ampere equation have 
appeared until recently. We make a distinction between numerical approaches with 
optimal transportation type boundary conditions Hj and the standard Dirichlet 
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boundary conditions (|D|). In the latter case, a number of numerical methods have 
recently been proposed for the solution of the Monge- Ampere equation. 

The early work by Benamou and Brenier [1] used a fluid mechanical approach 
to compute the solution to the optimal transportation problem. The work by 
Oliker and Prussner |39] presents a discretization that converges to the Aleksandrov 
solution in two dimensions. This is a very early method which solves a problem 
with only about a dozen grid points. 

For the Dirichlet problem treated herein, a series of papers have recently ap- 
peared by two groups of authors, Dean and Glowinski [HJ [T31 [M] and Feng and 
Neilan [TTJ HH] . The methods introduced by these authors perform best in the reg- 
ular case and can break down in the singular case. See [5] for a more complete 
discussion. 

We also mention the works j35j , in the periodic case, and [IS] for applications to 
mappings. The method of [46] treats the problem of periodic boundary conditions 
in an odd-dimensional space. 

1.4. Numerical challenges. When the conditions for regularity are satisfied, clas- 
sical solutions can be approximated successfully using a range of standard tech- 
niques (see, for example, works such as [H2 HH |2H and [TTJ, [H]). In an earlier 
work [5], we studied a simple finite difference discretization, which was accurate 
and fast on smooth solutions. However, for singular solutions, standard numeri- 
cal methods break down: either by becoming unstable, poorly conditioned, or by 
selecting the wrong (non-convex) solution. 

Weak solutions. For singular solutions, the appropriate notion of weak (viscosity 
or Aleksandrov) solution must be used. Numerical methods have been developed 
which capture weak solutions: Oliker and Prussner, in |39j . presented a method 
that converges to the Aleksandrov solution. The second author [37] introduced a 
wide stencil finite difference method which converges to the viscosity solution. Both 
of these methods were restricted to two dimensions. 

Convexity. The convexity constraint is necessary for both uniqueness and stability. 



In particular, the equation ( MA ) fails to be elliptic if u is non-convex (see subsec- 



tion 2.5), so instabilities can arise if the convexity condition (|Cj) is violated. Any 



approximation of ( MA ) requires some selection principle to choose the convex so- 
lution. This selection principle can be built into the discretization, as in [37] . or 
built into the solution method, as in [5]. 

Accuracy. The convergent monotone scheme of [37] uses a wide stencil and the 
accuracy of the scheme depends on the directional resolution, which in turn depends 
on the width of the stencil. As we demonstrate below, for highly singular solutions 



such as ( 13 ), the directional resolution error can dominate. However, it is unrealistic 
to expect high accuracy for singular solutions. Experimental results on singular 
solutions using standard finite differences, which is formally 0{h 2 ), produce results 
which are only accurate to O(vh); see [5|. 

Fast solvers. Previous work by the authors and a coauthor [5] investigated fast 



solvers for (MA I. An explicit method was presented which was moderately fast, 



independent of the solution time. For regular solutions, a faster (by an order of 



magnitude) semi-implicit solution method was introduced (see subsection 4.2 1, but 



this method was slower (by an order of magnitude) on singular solutions. 
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2. Analysis and weak solutions 
In this section we review regularity results and background analysis. 



Example solutions include (111 and (13 1, below. The first is a viscosity solu- 



tion, defined in subsection 2.2 The second is an Aleksandrov solution, defined in 
Isubscction 2.31 

The regularity conditions required for classical solutions are reviewed in the sec- 
tion below. The weaker notion of solution, viscosity solutions, allows the continuous 
function / to be zero. The even weaker geometric notion of Aleksandrov solutions 
allows / to be a nonnegative measure. 

2.1. Regularity. Regularity results for the Monge- Ampere equation have been 
established in (SH3HQZI]. We refer to the book [55] for the following result. 

Under the following conditions, the Monge-Ampere equation is guaranteed to 
have a unique C 2 ' a solution. 

{The domain, fl, is strictly convex with boundary 3D, € C 2 ' a . 
The boundary values, g € C 2a {dVl). 
The function, / € C a (fl) is strictly positive. 



Remark 1. In the extreme case, with f(x) = for all x € f2, the equation ( MA ) , (|C|) 
reduces to the computation of the convex envelope of the boundary conditions [361 
|3"8] . In this case, solutions may not even be continuous up to the boundary and 
can also be non-differentiable in the interior. 

Remark 2. While it is usual to perform numerical solutions on a rectangle, regu- 
larity can break down in particular convex polygons [451 141) . 

2.2. Viscosity solutions. We recall the definition of a viscosity solution [12"] . 
which is defined for the Monge-Ampere equation in |25) . 

Definition 1. Let u € C{Q) be convex and / > be continuous. The function 
u is a viscosity subsolution (supersolution) of the Monge-Ampere equation in Q, if 
whenever convex (f> e C 2 (S!) and x e ft are such that (u — <j>)(x) < (>)(u — 4>){xq) 
for all x in a neighborhood of Xq, then we must have 

det(D 2 4,(x )) > (<)f(x ). 

The function u is a viscosity solution if it is both a viscosity subsolution and su- 
persolution. 



Example 1 (Viscosity solution). This example will be computed in sections M7 



Consider (MA I in two dimensions, with solution, u, and data, /, given by 



«(x) = ^((|x|-l)+) 2 , /(x)=(l-ji 

Here we use the notation h + = max(/i, 0). (In three dimensions, the function / is 



different; see section 7). The function u is a viscosity solution but not a classical 
C 2 solution of the Monge-Ampere equation. 

We verify the definition of a viscosity solution. This only needs to be done at 
points where |xo| = 1 (since u is locally C 2 away from this circle). We note that / 
is equal to zero on this circle. 

We begin by checking convex C 2 functions </> < u with </>(x ) = m(x ) = (that 
is, u — <p has a local minimum here). Since Vu(x ) = 0, we require V0(x o ) = 
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as well. Since u is constant in part of any neighborhood of Xo, any convex <p must 
also be constant in this part of the neighborhood in order to ensure that u — 4> 
has a local minimum. This means that <fi has zero curvature in some directions so 
that detZ) 2 </>(xo) = 0, as required by the definition of the viscosity solution. We 
conclude that u is a supersolution of the Monge- Ampere equation. 

We also need to check functions <f> > u with </>(xo) = u(xq) = (so that u — <j> 
has a local maximum) . Since <f> is convex, it will automatically satisfy the condition 
det D 2 (/)(xq) > and we conclude that u is a subsolution. 

2.3. Aleksandrov solutions. Next we turn our attention to the Aleksandrov so- 
lution, which is a more general weak solution than the viscosity solution. Here 
/ is generally a measure [25]. We begin by recalling the definition of the normal 
mapping or subdifferential of a function. 

Definition 2. The normal mapping (subdifferential) of a function u is the set- 
valued function du defined by 

du(x ) = {p | u(x) > u(xo) + p ■ (x — Xo), for all x € £1} . 

For a set V C fi, we define du(V) = [J du(x). 

xev 

Now we want to look at a measure generated by the Monge- Ampere operator. 

Definition 3. Given a function u £ C(Q), the Monge- Ampere measure associated 
with u is defined as 

»(y) = \du{v)\ 

for any set V C fl, where \B\ is the measure of the set B. 

This measure naturally leads to the notion of the generalized or Aleksandrov 
solution of the Monge- Ampere equation. 

Definition 4. Let /i bea Borel measure defined in a convex set O e t d . The 
convex function u is an Aleksandrov solution of the Monge- Ampere equation 

det(L> 2 M) = fi 

if the Monge- Ampere measure associated with u is equal to the given measure /x. 

Example 2 (Aleksandrov solution). As an example, we consider the cone and the 
the scaled Dirac measure 



u(x) = |x|, n(V) =tt / <J(x)dx. 

Jv 

We verify from the definition that u, [i is an Aleksandrov solution of the Monge- 

Ampere equation. (Since [i is a measure, we cannot interpret a as a viscosity 

solution of Monge- Ampere.) The subdifferential du is given by 

du(x) - 

where Bi = {xEf 1 | |x| < l}. Then the associated Monge- Ampere measure will 
be 

\du(V)\ = r ° o € ^=nJ v S(x)dx = ^V). 
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2.4. A PDE for convexity. The convexity constraint ([C]) is necessary for unique- 
ness. For example, in two dimensions, ignoring boundary conditions, in the absence 



of the convexity constraint, —u will also be a solution of (MA I whenever u is a so- 
lution. 

For a twice continuously differentiable function u, the convexity restriction ([C]) 
is equivalent to requiring that the Hessian, D 2 u, be positive definite. Since we wish 
to work with less regular solutions, (JO) can be enforced by the equation 

Xi{D 2 u) >0, 

understood in the viscosity sense [36, 38 , where \i(D 2 u) is the smallest eigenvalue 
of the Hessian of u. 

The convexity constraint can be absorbed into the operator by defining 

d 

(3) det+(M) = n A i" 

where M is a symmetric matrix, with eigenvalues, Ai <...,< A„ and 

x + = max(a;,0). 



Using this notation, (MA|,(0 becomes 



(MA)+ det + (D 2 u(x)) = f(x), for x in il. 

Remark 3. Notice that there is a trade-off in defining pi: the constraint (lei) is 
eliminated but the operator becomes non-differentiable near singular matrices. 

2.5. Linearization and ellipticity. The linearization of the determinant is given 

by 

V det(M) • N = trace {M ado N) 

where M a( g is the adjugate [42], which is the transpose of the cofactor matrix. The 
adjugate matrix is positive definite if and only if M is positive definite. When the 
matrix M is invertible, the adjugate M a( u satisfies 

(4) M adj = det(M)M-\ 

We now apply these considerations to the linearization of the Monge-Ampere 
operator. When u € C 2 we can linearize this operator as 

(5) V M det(L> 2 u) ■ v = trace ((D 2 u) adj D 2 v) . 
Example 3. In two dimensions we obtain 

Vm det(D u)v = u xx v yy + u yy v xx - 2u xy v xy , 

which is homogeneous of order one in D 2 u. In dimension d > 2, the linearization 
is homogeneous order d — 1 in D 2 u. 

The linear operator 

L[u] = tra,ce(A(x)D 2 u) 
is elliptic if the coefficient matrix A(x) is positive definite. 

Lemma 1. Let u € C 2 . The linearization of the Monge-Ampere operator, (OH), is 
elliptic if D 2 u is positive definite or, equivalently, if u is (strictly) convex. 

Remark 4. When the function u fails to be strictly convex, the linearization can be 
degenerate elliptic, which affects the conditioning of the linear system fcfy. When 
the function u is non-convex, the linear system can be unstable. 
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The definition of a nonlinear elliptic PDE operator, which follows, generalizes 
the definition of a linear elliptic operator. In addition, it allows for the operator to 
be non-differentiable. 

Definition 5. Let F{M) be a continuous function defined on symmetric matrices, 
which we regards as a PDE operator. Then F(M) is elliptic if it satisfies the 
monotonicity condition 

F(M) < F{N) whenever M < N, 

where for symmetric matrices M < N means x T Mx < x T ' Nx for all x. 

Example 4. The operator det (M) is a non-decreasing function of the eigenvalues, 
so it is elliptic. 

3. Convergent discretization of the Monge-Ampere equation 

In this section, we produce a discretization of the Monge-Ampere equation in 
two and higher dimensions, and prove that solutions of the scheme converge to the 
viscosity solution. This discretization is different from the one in |37j . which was 
restricted to two dimensions. 

We also prove convergence of Newton's method to the solution of the scheme. To 
this end, we modify the original monotone scheme from a non-differentiable scheme 
to a regularized (but still monotone) scheme. 

3.1. The standard finite difference discretization. We begin by discussing the 
standard finite difference discretization of the Monge-Ampere equation. This is ob- 
tained by simply discretizing each of the operators using standard finite differences 
as in, for example, [5|. 

Since this discretization alone does not enforce the convexity condition ([C]), it 
can lead to instabilities. In particular, Newton's method can become unstable if 
this discretization is used. 

There is no reason to assume that the standard discretization converges. In 
fact, the two dimensional scheme has multiple solutions. In 5j this discretization 
was used, but the solvers were designed to select the convex solution. In addition, 
the solution methods used in [5] do not generalize to three and higher dimensions: 
enforcing convexity of the solution is much more difficult when there are more that 
two possible eigenvalues of the Hessian. 

3.2. Eigenvalues of the Hessian in two dimensions. In two dimensions, the 
largest and smallest eigenvalues of a symmetric matrix can be represented by the 
Rayleigh-Ritz formula 

Ai (A) = min v T Av, A? (A) = tok^v t Av. 

{ ' 1-1=1 V ' 1-1=1 

This formula was used in |37] to build a monotone scheme for the Monge-Ampere 
operator, which is the product of the eigenvalues of the Hessian. 

In higher dimensions, the formula above does not generalize naturally. Instead, 
we use another characterization, which applies to positive definite matrices. 
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3.3. A variational characterization of the determinant. In this section we 
establish a matrix analysis result that is used to build a monotone discretization 
of the Monge- Ampere operator. Our result below is closely related to Hadamard's 
inequality, which states that, for a symmetric positive definite matrix A, 



dct(A) < Yl 



i=l 

with equality when A is diagonal. 

Consider an arbitrary symmetric positive definite matrix A. We can characterize 
the determinant of A as follows. 

Lemma 2 (Variational characterization of the determinant). Let A be a d x d 
symmetric positive definite matrix with eigenvalues Xj and let V be the set of all 
orthonormal bases of R d : 

V = {{v u ...,v d )\v i € M, d , Vi J_ v 3 if i ^ j, \\ Vj \\ 2 = 1}. 

Then the determinant of A is equivalent to 

d d 

1 [ A, = min I [ vj Au-j. 

Proof. Since A is symmetric and positive definite, we can find a set of d orthonormal 
eigenvectors Vj. 

Any (yi, . . . , Vd) € V, can be expressed as a linear combination of the eigenvec- 
tors: 

d d 

Vj = ^CjkVk = ^2(vjv k )v k . 
fe=l fe=l 

Since the Vj and Vj are both orthonormal, we can make some claims about the 
coefficients Cjk- 



d / d \ / d \ 

Y c % = Y c Jk v I Y c J m 
fe=i \fe=i / \i=i / 



= v J v i = !. 



d I d \ 

Y c ?fc = v k \Y v j v J Vk = v ^ vk = L 

We can use these results to compute 

d d d / d \ 

log J] ^ = £ log (i/J^-) = 2 log ^ s 2 fc A fe 

.7=1 j=l J=l u=l / 

Using Jensen's inequality, we conclude that 

rf d d 

log JJ if^ > £ E s 2 ^ lo s( A fc) 



j=l .7=1 fc=l 

d / d \ d 

=y \Y c % iog(A fe )-io g n^. 

fc=l \j=l / J=l 
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Since the logarithmic function is increasing, we conclude that 

d d 

IKM>n A ; 

3=1 J = l 

with equality if the Va are identical to the eigenvectors Vj . Thus we conclude 

d d 

Y\ X, = min \\vJAua. D 

J = l 3=1 

We use Lemma [2J to characterize the determinant of the Hessian of a convex C 2 
function <f>. We can express this in terms of second directional derivatives of <j> as 
follows: 

d d „2 1 

det(D 2 (b)= min \\v^D 2 ^v j = min TT —%. 
(^,...«)ev lj ; 3 {v u ...,v d )av i - i -dv : j 2 

3=1 3=1 J 



The convexified Monge- Ampere operator (M A) + is represented by 



,; / a 2 -' vl 



det+( J D 2 0) = min 17 ^- 



2 I 
J=l \""J/ 

Remark 5. This characterization of the Monge- Ampere operator remains valid even 
when <fi is not strictly convex, in which case the value is zero. 

3.4. Wide stencil schemes. Wide stencil schemes are needed to build consistent, 
monotone discretizations of degenerate second order PDEs [7J [33J 131] ■ Wide stencil 
schemes were built for the two dimensional Monge- Ampere equation in [37] . A wide 
stencil discretization of the convex envelope was given in |36j . 

To discretize the Monge- Ampere operator on a finite difference grid, we approx- 
imate the second derivatives by centered finite differences; this is the spatial dis- 
cretization, with parameter h. In addition, we consider a finite number of possible 
directions v that lie on the grid; this is the directional discretization, with param- 
eter d9. We denote this set of orthogonal vectors by Q . Then we can discretize the 
convexified Monge- Ampere operator as 

d 

(MA) M MA M [u]= min TT (V v . v .u) + 

{vi...v d }eG ~ 
j=i 

where T> vv is the finite difference operator for the second directional derivative in 
the direction u, which lies on the finite difference grid. These directional derivatives 
are discretized by simply using finite differences on the grid: 

T^wUi — — 2 — ( u ( x i + vh) + u(xi — vh) — 2u{xi)) . 

\v\ h 2 

Depending on the direction of the vector u, this may involve a wide stencil. At 
points near the boundary of the domain, some values required by the wide stencil 



will not be available (Figure 1 ). In these cases, we use interpolation at the boundary 
to construct a (lower accuracy) stencil for the second directional derivative; see [37] 
for more details. 

Since the discretization considers only a finite number of directions v 1 there will 
be an additional term in the consistency error coming from the directional resolution 
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Figure 1. Wide stencils on a two-dimensional grid. 



dO of the stencil. This angular resolution will decrease and approach zero as the 
stencil width is increased. In practice, we use fixed relatively narrow stencils for 
most computations. However, for singular solutions the directional resolution error 
can dominate. 

We also note that the discretization we have just described is genuinely distinct 
from the two-dimensional wide stencil discretization described in [37] . For example, 
we can consider the function u(x,y) — x 2 + y 2 + x 2 y 2 and discretize the Monge- 
Ampere operator using a 9-point stencil. This allows us to choose from the set of 
directions 

{(1,0), (0,1), (1,1), (1,-1)}. 

Using the two-dimensional characterization of the Monge-Ampere equation (re- 
called in subsection 3.2), the monotone discretization produces 



min T>,, 



,u 



( 



max T>,, 



, u 



V\ J \ VI 

On the other hand, our new discretization has the value 

min {T> VlVl uV V2U2 u} = 4. 

Vl -L V2 



2(2 + h 2 



3.5. Regularized monotone discretization. The monotone discretization we 
have just described may not be differentiable at points where the minimum is 
attained along more than one direction z/, or at points where the value is zero. 
Since we need to differentiate the operator when we build fast solvers using Newton's 
method in |scction 5[ we wish to regularize this discretization. To ensure convergence 
to the viscosity solution, we need to make the regularization monot one. 

One way to do this is to notice that the non-differentiability of (MA) M arises 
only from the operations of max and min. Thus if we regularize each of these 
operations in a monotone way, we can reconstruct a regularized version of (MA) 
by substitution. 
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With that in mind, we define 
<5 



max 

• 6 

mm 



(a,b) = \ (a + b + y/(a - 6) 2 + <5 2 ) 
(a 7 b) = ^(a + b-V(a-by + 5^ 



Clearly max — » max and min — > min as 5 — > 0. Moreover, these functions are 
differentiable and non-decreasing in each variable. 

Now we can build up the regularized operator as follows. Define 

{V ul ,u) + ~ = max* (T> vv u, 0) 



and replace each term in the product in {MA) with its regularized version as 



above. Next, the minimum of the product in {MA) is regarded as a succession 
of minimums, each of which is replaced by its regularized version. 
The resulting discretization is denoted 

(MA) S MA s [u], 

It is a smooth function of u,, strictly increasing in each of the T> v k v kUi, and converges 



to the original discretization {MA) as 5 — > 



3.6. Convergence of finite difference schemes. In order to prove convergence 
of the solutions of our finite difference schemes to the unique viscosity solution 



of ( MA I , we will rely on a framework developed by Barles and Souganidis [3] and 
extended in [33] . 

The framework of [3] gives easily verified conditions for when approximation 
schemes converge to the unique viscosity solution of a PDE. 

Theorem 3 (Convergence of Approximation Schemes). Consider a degenerate el- 
liptic equation, for which there exist unique viscosity solutions. A consistent, stable 
approximation scheme converges uniformly on compact subsets to the viscosity so- 
lution, provided it is monotone. 

While the previous theorem gives conditions for convergence, it does not pro- 
vide a method for building monotone schemes or verifying when the schemes are 
monotone. In |35| . a framework for building and verifying the monotonicity of finite 
difference schemes was established. 

We recall that a finite difference equation has the form 

F l [u] =F i {u h u i -u j \ jl ti). 

This allows us to characterize a degenerate elliptic (monotone) scheme as fol- 
lows [55] , 

Definition 6. The scheme F is degenerate elliptic if F l is non-decreasing in each 
variable. 

We recall Theorem 3 from |35j , which provides a simple way of verifying both 
monotonicity and stability. 

Theorem 4 (Verifying Monotonicity and Stability). A scheme is monotone and 
non-expansive in the t°° norm if and only if it is degenerate elliptic. 

The notion of degenerate ellipticity is also enough to guarantee the existence of 
a unique solution to a scheme, as proved in Theorem 8 of [33] , 
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Theorem 5 (Existence and Uniqueness of Solutions for elliptic schemes). A proper, 
locally Lipschitz continuous degenerate elliptic scheme has a unique solution which 
is stable in the £°° norm. 



Given these general results, the work in proving that a discretization of (MA) 



converges is reduced to verifying two conditions: consistency and degenerate cllip- 
ticity. 

Remark 6 (Convergence rates). While the formal accuracy of the scheme can be 
determined by Taylor series applied to smooth test functions, the convergence the- 
orem guarantees only uniform convergence. This is to be expected since, in general, 
viscosity solutions can be singular. The power of the convergence result is that it 
applies even in the singular case. In general, the observed accuracy depends on 
both the regularity of the solution and the choice of discretization, with observed 
values going from 0(h 2 ) (for C 4 solutions) to 0(Vh); see [5]. 

3.7. Proof of convergence. 



Theorem 6 (Convergence to Viscosity Soluti on) 



viscosity solution. Solutions of the schemes (MA) M 
cosity solution of (MA I as h,d9,S — > 0. 



(MAy 



Let t he PDE (MA) have a unique 
converge to the vis- 



Proof. The convergence follows from verifying consistency and degenerate ellipticity 
as in 37]. This is accomplished in Lemmas 7|8 □ 



Stability in £°° for solutions of the schemes follows from the fact that the schemes 
are degenerate elliptic and the results of Theorem [5] 



Le mma 7 ( Degenerate Ellipticity). The finite difference schemes given by (MA) M 
and 



(MA) are degenerate elliptic. 



Proof. From their definition, the discrete second directional derivatives T> vv u are 
non-decreasing functions of the differences between neighboring values and reference 
values, Uj — Uj, where Uj is one of the neighboring values in the direction v. The 



scheme (MA) is a non-decreasing combination of the operators min and max 



applied to the degenerate elliptic terms T> vv u 1 so it is also degenerate elliptic. 



We recall from th e construction of the scheme in subsection 3.5 that the regular- 
ized scheme (MA) comes from replacing the operations of min and max in (MA) 
by a non-decreasing regularization of these operations. Thus the combined scheme 
is also degenerate elliptic. □ 



We also require the schemes (MA) 



and 



Monge- Ampere equation. We prove consistency of (MA) since consistency of (MA) 
is a special case. 



(MA) M to be consistent wi th the 



\M 



Definition 7. The scheme MA • ■ is consistent with the equation (MA) at xq if 
for every twice continuously differentiable function <j)(x) defined in a neighborhood 

of Xq, 



MA 



h,dd,S 



[4>}(x ) -> MA[4>](x ) as h,d6,6 -*■ 0. 



The global scheme defined on f2 is consistent if this limit holds uniformly for all 
x e Q. 
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Lemma 8. Let xq € fl be a reference point on the grid and 4>(x) be a twice 
continuously differentiable function that is define d and co nvex in a neighborhood 
of the grid. Then the scheme M J 4 <5 [0] defined in (MA) S approximates the PDE 
MA[(j)} with accuracy 

MA 5 [cj)}(x ) = MA[cf>}(x ) + 0(h 2 + d8 + r(6)) 

where the function r(S) converges to zero as S — > 0. 

Proof. From a simple Taylor series computation we have 

V vu (t){x ) = <t> vv (x ) +0(h 2 ). 

We also recall that in |subscction 3.5| we regularized the second directional derivatives 
to obtain 

V s vu <f>(x ) = mfflc{D w ^o), 0} + O(S) = <j> vv (x a ) + 0(h 2 + S). 

Here the maximum has no effect since we are considering convex <p. 
Now we recall that the Monge- Ampere operator can be expressed as 



d 
veV 






where the Vj are orthogonal unit vectors, which may not be in the set of grid vectors 
Q. We can then choose a set of vectors 

v + dv 

6 6, 



I v + dv I 



so that each remainder \dvj\ — O(d0). 

Now we consider the discretized problem 



mmT[V s cj>(x ) = mmT\V VjUj <t>{x Q )+0{5) 

- II ^(Vj+dv^ivj+dvj^ixo) + 0(S) 






j=l 



y-r (^ + dv j ) T D 2 (j)(x )(v :j + (%) ^ 2 

M h+^-i 2 r ' ' 



fj vJD 2 ^(x )v J + 0(h 2 +d6 + S) 

3=1 

d 

min FT 0„ „ (z ) + 0(^ 2 + dfl + S). 



J' = l 
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In addition, since the set of grid vectors Q is a subset of the set of all orthogonal 
vectors V, we find that 

d d 

mmTXvl v Mx Q ) > mmT[V VjVj (P(xo) + 0(S) 

3 = 1 J=l 

d 

= min TT <j> vv (x ) + 0(h 2 + S). 
We conclude that 

d d 

min T\ Vl v .(j){x Q ) = min T\ <j> vu (x ) + 0{h 2 + d9 + 6), 

which is precisely the characterization of the Monge- Ampere operator given in |Thc-| 
lorem 2 1 

Finally, we replace the minimum in the above scheme with a smooth function, 



which converges uniformly as 8 — > by construction (subsection 5.1). Thus the 
resulting scheme will satisfy 

MA & [<f>]{xo) = det{D 2 4){x Q )) + 0(h 2 + d6 + r(S)). □ 

4. A SEMI-IMPLICIT SOLUTION METHOD 



Any discretization of ( MA I leads to a system of nonlinear equations that must 
be solved in order to obtain the approximate solution. Newton's method requires 
a good initial value to converge. 



In this section, we describe a semi-implicit solution method for (MA). One 
iteration of this method will be used to compute the initial value for Newton's 
method. 

First we describe the fully explicit method. 

4.1. Explicit solution methods for monotone schemes. Using a monotone 
discretization F[u] of the Monge- Ampere operator, the fixed point iteration, 

u n+l = u n + di ( f [ u r,] _ ^ 

is a contraction in £°° , provided dt is small enough. This iteration corresponds to 
solving the parabolic version of the equation using a forward Euler discretization. 
Explicit iterative methods have the advantage that they are simple to implement, 
but the number of iterations required suffers from the well known CFL condition 
(which applies in a nonlinear form to monotone discretizations, as explained in [35]). 
This approach is slow because for stability it requires a small time step dt, which 
decreases with the spatial resolution h. The time step, which can by computed 
explicitly, is 0(h 2 ). This was the approach used in [57] . 

4.2. A semi-implicit solution method. The next method we discuss is a semi- 
implicit method, which involves solving a Poisson equation at each iteration. 

In [5] we used the identity Q to build a semi-implicit solution method. We 
showed that the method is a contraction, but the strictness of the contraction 
requires strict positivity of /. In practice, this meant that the iteration was fast 
for regular solutions, but became slower than the explicit method when / was zero 
in large parts of the domain. 
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We begin with the following identity for the Laplacian in two dimensions, 



|Au| = ^(Au) 2 = ^ju 2 xx + u 2 yy + 2u xx u yy . 
If u solves the Monge- Ampere equation, then 



| A«| = Ju 2 xx + u 2 + 2u 2 + 2/ = V \D 2 u\ 2 + 2f. 



'.iy 



This leads to a semi-implicit scheme for solving the Monge- Ampere equation, used 
in 0. 



(6) Au n+1 = y2f+\D 2 u n \ 2 . 

To generalize this to K d , we can write the Laplacian in terms of the eigenvalues of 
the Hessian: Am = J2i=i Xi[D 2 u]. Taking the <i-th power and expanding gives the 
sum of all possible products of d eigenvalues 

d 

(Au) d = d\l[X i + P(X 1 ,...,X d ), 

where -P(A) is a d-homogeneous polynomial, which we will not need explicitly. 
The result is the semi-implicit scheme 

(7) Au n+1 = (d\f + P(X t [D 2 u n }, . . . , X d [D 2 u n ])) 1/d . 
A natural initial value for the iteration is given by the solution of 

(8) Au° = (d!/) 1/d . 

5. Newton's method 



We perform the damped Newton iteration 
to solve the discretized equation 



u n+l = u n _ av r. 



MA M [u] = f, 

where the damping parameter a, < a < 1, is chosen at each step to ensure that 
the residual \\MA M {u n ) — f\\ is decreasing. 
The corrector v n solves the linear system 

(9) (\7 U MA M [u n ]) v n = MA M [u n ] - f. 

Here the left hand side is our notation for the Jacobian matrix of the scheme. The 
Jacobian matrix for the monotone discretization is obtained by using Danskin's 
Theorem [B] and the product rule. 

V u MA M [u] = J2 diag ( J] -Du M u\ V„. v . 



where the v* are the directions active in the minimum in 



(MA) M 
In order for the linear equation |9]) to be well-posed, we require the coefficient 
matrix to be positive definite. As observed in Lemma [T] this condition can fail if 
the iterate u n is not strictly convex. 
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5.1. Convergence of Newton' s Method for the regularized discretization. 

While the discretization (MA) M is designed to be robust enough to converge to sin- 



gular solutions of (MA), we wish to compute these solutions quickly using Newton's 



Method. Newton's Method requires that the objective function be differentiable (or 
at least dir ectionally differentiable) with an invertible Jacobian matrix. However, 
the scheme (MA) M is non-differentiable and, even at points of differentiability, 
it can have a degenerate Jacobian matrix. The latter possibility corresponds to 
non-strictly convex solutions, which have degenerate Hessians. Newton's method 
applied to this scheme can break down for singular solutions. 



The objective of this section is to show that the regularized scheme {MA) 



is 



designed so that even singular solutions can be computed using Newton's method. 

Theorem 9 (Newton's Method for the Discretized Monge- Ampere Equation). Sup- 
pose the PDE ( MA I has a unique v iscosity solution. Then Newton's method for the 
discretized system given by (MA) converges quadratically. 



In order to prove this result, we recall a standard result on the convergence of 
Newton's method for a system of equations [31]. 

Lemma 10 (Newton's Method for a System of Equations). Consider a system of 



equations F[u] — where the operator F : K - 
Suppose the following conditions hold: 

(1) A solution u* e U exists. 

(2) VF : U -> R NxN is Lipschitz continuous. 

(3) VF(u*) is non-singular. 

Then the Newton iteration 



u n+1 = u n - VF{u n )- l F(u n ) 
converges quadratically to u* if u° € U is sufficiently close to u* 



and let U C R be open. 



Remark 7. In the proof below we use the fact that the discretization (MA)"\ is 
degenerate elliptic, which leads to a positive definite and nonsingular Jacobian VF 
in the Newton iteration. 

Proof of Theorem [9[ For any fixed grid, the discretized system of equations has a 
solution, as es tablished in Theorem [6] 

The scheme (Ad A) 5 is smooth in u and the Jacobian is therefore locally Lipschitz 
continuous. 

By construction, the discrete Monge-Ampere operator is strictly increasing in 



each of the discrete second directional derivatives (subsection 3.5 1. Thus the Jaco- 
bian will have the form 



V u MA s [u] = J2 X^ifcM^fc 

v k eg j=i 

where each of the Ajk(u) is a positive definite diagonal matrix. The Jacobian is 
negative definit e an d thus invertible. 
By Theorem [lO 



Newton's method converges for the discretized system (MA) 



□ 
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5.2. Initialization of Newton's method. Newton's method requires a good ini- 
tialization for the iteration. Since we need the resulting linear system to be well- 
posed, it is essential that the initial iterate: (i) be convex, (ii) respect the boundary 
conditions, (hi) be close to the solution. In addition, the computational cost of 
initialization should be low. 

In order to achieve these conditions we use one step of the semi-implicit scheme ufh 
to obtain a close initial value. This amounts to solving ^ along with consistent 
Dirichlet boundary conditions (JDJ). As necessary, we then convexify the result using 
the method of |36 (this last step was only needed for the most singular solution). 
Since both these steps can be performed on a very coarse grid and interpolated 
onto the finer grid, the cost of initialization is low. 

5.3. Preconditioning. In degenerate examples, the PDE for v n Q may be de- 
generate, which can lead to an ill-conditioned or singular Jacobian. To get around 
this problem, we regularize the Jacobian to make sure the linear operator is strictly 
negative definite; this will not change the fixed points of Newton's method. We 
accomplish this by replacing the second directional derivatives u vv with 



u vv = max{« w ,e}. 

Here e is a small parameter. In the computations of 
1(T 8 . 



section 6 



we take e = ^ 



6. Computational results in two dimensions. 

In this section, we summarize the results of a number of two-dimensional ex- 
amples computed using the methods described in this paper. In particular, we are 
interested in comparing the computation time for Newton's method with the time 
required by the methods proposed in [5]. 

We perform the computations using a 17 point stencil on an N x N grid on a 
square. 

We solved all linear systems using the MATLAB backslash operator, which per- 
forms an LU decomposition of the sparse systems arising in Newton's method. 
Naturally, computations could be made faster by using a compiled programming 
language or a more efficient linear solver. However, this implementation is sufficient 
to demonstrate the efficiency of our method. 

6.1. Four representative examples. We have tested the monotone scheme on a 
number of examples of varying regularity. To illustrate these, we present detailed 
results for four representative examples. 

To define these exact solutions, we first write x = (x, y) and xo = (.5, .5) for the 
center of the domain. 

The first example, which is smooth and radial, is given by 

(10) u(x) = exp (^pj , /(x) = (l + |x| 2 ) exp (|x| 2 ) . 
The second example, which is C , is given by 

(11) u(x) = i((|x-x |-0.2)+) 2 , /( x )=( 1 "^ Co | 
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The third example is twice diffcrentiable in the interior of the domain, but has 
an unbounded gradient near the boundary point (1, 1). The solution is given by 

(12) u(x) = - A /2-|x| 2 , /(xH2(2-|x| 2 )~ 2 . 

The final example is the cone, which was discussed in |subsection 2.3| It is 
Lipschitz continuous. 

(13) u(x) = v/|x-x Q |, / = M = tt c5 Xo . 

In order to approximate the solution on a grid with spatial resolution h using 
viscosity solutions, we approximate the measure [i by its average over the ball of 
radius h/2, which gives 



f h = 



A/ti 2 for |x-x | < h/2, 
otherwise. 



6.2. Visualization of solutions and gradient maps. In Figure 2 the solutions 
and the gradient maps for the first three representative examples are presented. 



For example (13), the gradient map is too singular to illustrate. To visualize the 



maps, the image of a Cartesian mesh under the mapping 

x \ / T> x u 

V ) V V v u 

is shown, where (T> x u, 2? y u) is the numerical gradient of the solution of the Monge- 
Ampere equation. The image of a circle is plotted for visualization purposes; the 
equation is actually solved on a square. For reference, the identity mapping is also 
displayed. 

In each case, the maps agree with the maps obtained using the gradient of the 
exact solution. 







Regularity of Solution 




Method 


C 2a 


C l,a 




C ' 1 (Lipschitz) 


Gauss-Scidcl 


Moderate 


Moderate 




Moderate 




(~ 0{M 1A )) 


(~ OiM 1 - 9 )) 




(~ G(M 2 )) 


Poisson 


Fast 


Slow-Fast 




Slow 




(~ 0(M 1A ) 


~ 0(M x - 4 )-blow-up) 


{- 


- C(M 2 )-blow-up) 


Newton 


Fast 


Fast 




Fast 




(~ (DIM 1 - 3 )) 


(~ OiM 1 - 3 )) 




(~ 0(M 13 )) 



Table 1 . Approximate computation time required by the Gauss- 
Seidel, Poisson, and Newton's methods for two-dimensional prob- 
lems of varying regularity. Here M = N 2 is the total number of 
grid points. 



6.3. Computation time. The computation times for the four representative ex- 



amples are presented in Table 2 The computations time are compared to those for 



the Gauss-Seidel and Poisson iterations described in |S]. In each case, the Newton 
solver is faster in terms of absolute solution time. 

|Tablc l] presents order of magnitude solution times. The order of magnitude 
solution time for Newton's method is independent of the regularity of the solution 
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-1 -1 



(a) 




(b) 





(d) 



0.04 



0.02- 





(f) 





(h) 



Figure 2. Solutio ns and m appings for the |(a)|(b)| identity map, 
(c)|(d) C 2 example, (c)|(f) C 1 example, and | (g) | (h)j example with 
blow-up. 
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and faster than both of the other methods. In particular, the computation time 
required by our Newton solver is roughly 0(M 13 ). This is really dependent on two 
different things: the number of Newton iterations required, which can increase as 
the grid is refined, and the speed of the linear solver. The use of a more efficient 
linear solver might further reduce the order of magnitude solution time for this 
method. 







C 2 Example ( 10 1 




N 


Newton 


CPU Time ( 


seconds) 




Iterations 


Newton 


Poisson 


Gauss-Seidel 


31 


3 


0.2 


0.7 


2.2 


G3 


6 


0.9 


1.9 


15.0 


127 


7 


4.9 


9.6 


236.7 


255 


7 


43.5 


52.6 


— 


361 


7 


172.8 


162.6 


— 
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N 


Newton 


CPU Time ( 


seconds) 




Iterations 


Newton 


Poisson 


Gauss-Seidel 


31 


4 


0.3 


1.1 


0.8 


03 


7 


0.8 


20.5 


9.5 


127 


11 


6.6 


256.8 


145.5 


255 


16 


61.5 


— 


— 


361 


20 


350.4 


— 


— 



Example with blow-up ([12]) 
N Newton CPU Time (seconds) 

Iterations Newton Poisson Gauss-Seidel 



31 


4 


0.3 


0.5 


0.8 


63 


4 


0.5 


2.9 


19.4 


127 


5 


3.7 


17.7 


293.3 


255 


7 


41.4 


128.2 


— 


361 


9 


184.1 


374.5 


— 



C - 1 (Lipschitz) Example ([l3l 
N Newton CPU Time (seconds) 

Iterations Newton Poisson Gauss-Seidel 



31 


9 


0.5 


5.3 


0.8 


63 


15 


1.4 


91.9 


21.5 


127 


32 


14.1 


1758.2 


373.9 


255 


34 


101.7 


— 


— 


361 


29 


280.2 


— 


— 



Table 2. Computation times for the Newton, Poisson, and Gauss- 
Seidel methods for four representative examples. 
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6.4. Accuracy. Finally, we present accuracy results for the four representative 
examples; see |Table 3[ We perform the computations using the monotone scheme 
on 9, 17, and 33 point stencils. The accuracy of the scheme is determined by a 
combination of the directional resolution (d0) error and the spatial discretization 
error. Widening the stencil, which has the effect of decreasing dO, improves the 
accuracy, as does increasing the number of grid points. We also compared the 
accuracy to standard finite differences using the results of [5] . (Solution times were 
longer using these methods). Standard finite differences are formally more accurate 
since there is no d6 error. The numerical results show that for the first two, more 
regular examples, the standard finite differences are more accurate. However, for 
the two more singular examples, the monotone finite differences are slightly more 
accurate. 



7. Computational results in three dimensions 

Next, we perform computations to test the speed and accuracy of Newton's 
method for three-dimensional problems. These computations are done using a 19 
point stencil on an iV 3 grid on the square [0, l] 3 . 

The methods of [5j were restricted to the two-dimensional Monge- Ampere equa- 
tion, so no computations were available for comparison in three dimensions. 

As before, we performed computations on three representative exact solutions of 
varying regularity. To simplify the following expressions, we denote a vector in K 3 

by 

x= (x,y,z), 
and let 

x = (.5, .5,.5) 
be the center of the domain. The first example solution is the C 2 function 

(14) «(x) = exp faf\ , /(x) = (l + |x| 2 ) exp (^f-) . 

The second example solution is the C 1 function 

(15) U (x) = ((|x-x o |-0.2)+) 2 /2, 

fl-r^ + r^^, |x-x |>0.2, 

/(x) = < l x - x "l |x-x | 2 ' ' Ul ' 

1 otherwise. 

The third example is the surface of a ball, which is differentiable in the interior of 
the domain, but has an unbounded gradient at the boundary: 

2 tl , „/ ,2\- 5 / 2 



(16) U (x) = -^3-|x| z , /(x)=3(3-|x 

Computation times and accuracy results for the three-dimensional examples are 
presented in |Table 4) 
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N 


Max Error, C 2 
9 Point 17 Point 


Example (10) 
33 Point 




F.D. 


31 


17.9 x 1(T 4 


8.9 x 10~ 4 


7.0 x 10~ 4 


7.14 


x 10~ 5 


63 


16.2 x 1CT 4 


5.1 x 10~ 4 


3.1 x 10~ 4 


1.73 


x 10~ 5 


127 


15.9 x 1CT 4 


4.6 x 10~ 4 


1.8 x 10~ 4 


4.3 


x 10~ 6 


255 


15.9 x 1(T 4 


4.4 x 10~ 4 


1.5 x 10~ 4 


1.1 


x 10 -6 


361 


15.9 x 1CT 4 


4.4 x 10~ 4 


1.5 x 10~ 4 


0.5 


x 10- 6 


N 


Max Error, C 1 
9 Point 17 Point 


Example (11) 
33 Point 




F.D. 


31 


3.0 x 10~ 3 


1.7 x 10~ 3 


1.5 x 10~ 3 


2.6 


x 10~ 4 


G3 


2.5 x 1CT 3 


1.0 x 10~ 3 


0.6 x 10~ 3 


1.5 


x 10~ 4 


127 


2.3 x 1CT 3 


0.8 x 10~ 3 


0.3 x 10- 3 


0.6 


x 10~ 4 


255 


2.2 x 1CT 3 


0.7 x 10~ 3 


0.3 x 10~ 3 




— 


361 


2.2 x 1CT 3 


0.7 x 10~ 3 


0.3 x 10~ 3 






N 


Max Error, Example with blow-up ( 12 ) 

9 Point 17 Point 33 Point Standard F.D. 


31 


1.7 x 1CT 3 


1.7 x 10~ 3 


1.7 x 10~ 3 


17.15 


x 10 -3 


G3 


0.9 x 10~ 3 


0.6 x 10~ 3 


0.6 x 10~ 3 


12.53 


x 10~ 3 


127 


0.8 x 10~ 3 


0.3 x 10~ 3 


0.2 x 10~ 3 


9.00 


x 10 -3 


255 


0.8 x 10~ 3 


0.3 x 10~ 3 


0.2 x 10~ 3 


6.42 


x 10 -3 


361 


0.8 x 10~ 3 


0.3 x 10~ 3 


0.2 x 10~ 3 


5.41 


x 10~ 3 


Max Error, C ' 1 (Lipschitz) Example (13) 
N 9 Point 17 Point 33 Point Standard F.D. 


31 


12 x 10~ 3 


3 x 10~ 3 


3 x 10" 3 


10 


x 10" 3 


63 


11 x 10~ 3 


3 x 10~ 3 


2 x 10~ 3 


6 


x 10~ 3 


127 


11 x 10~ 3 


4 x 10~ 3 


2 x 10- 3 


3 


x 10~ 3 


255 


11 x 10~ 3 


4 x 10~ 3 


1 x 10~ 3 




— 


361 


11 x 10~ 3 


4 x 10~ 3 


1 x 10" 3 




— 



Table 3. Accuracy of the monotone scheme for different stencil 
widths and for standard finite differences on four representative 
examples. 



8. Conclusions 

A fast, convergent finite difference solver for the elliptic Monge- Ampere equation 
was built, analyzed, and implemented. Computational results were presented using 
two- and three-dimensional exact solutions of varying regularity, from smooth to 
non-differentiable. 

A monotone discretization was built using a method which applied arbitrary 
dimensions. A proof of convergence of the finite difference approximation to the 
unique viscosity solution of the equation was given. 

The discretized equations were solved using Newton's method, which is fast, 
experimentally 0(M 13 ) where M is the number of data points, independent of 
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C 2 


Example ( 14 1 


N 


Max Error 


Iterations 


CPU Time (s) 


7 


0.0151 


2 


0.1 


11 


0.0140 


3 


0.1 


15 


0.0132 


5 


0.5 


21 


0.0127 


6 


3.6 


31 


0.0125 


5 


34.7 



C 1 Example ([l5l 



N Max Error Iterations CPU Time (s) 



7 


0.0034 


1 


0.02 


11 


0.0022 


1 


0.06 


15 


0.0019 


1 


0.17 


21 


0.0020 


2 


1.42 


31 


0.0019 


2 


16.70 



Example with Blow-up ( 16 1 



N Max Error Iterations CPU Time (s) 



7 


9.6 x 10" 3 




1 


0.02 


11 


5.3 x 10^3 




3 


0.12 


15 


4.7 x 10^3 




3 


0.34 


21 


4.3 x 10~ 3 




6 


3.65 


31 


3.9 x 10~ 3 




8 


54.70 


Table 4. Maximum error 


and 


computation 


time on three repre 


sentative three-dimensional 


exa 


mples. 





the regularity of the solution. The implementation of Newton's method was signifi- 
cantly (orders of magnitude) faster than the two other methods used for comparison. 
A close initial approximate solution was provided using a low cost method based on 
performing one step of a previously established solution method, then taking the 
convex envelope as needed. A proof of convergence of Newton's method was also 
provided. 

The solver presented here used a novel discretization in general dimensions, ac- 
companied by a fast solution method. In terms of solution time and stability, the 
resulting solver is a significant improvement over existing methods for the solution 
of the elliptic Monge- Ampere equation. This conclusion is valid for both regular 
and singular solutions. On smooth solutions the accuracy was lower than when us- 
ing standard finite differences, but on singular solutions the accuracy was slightly 
better. 

After the present paper was submitted, we succeeded in improving the accuracy 
of the discretization by applying the monotone method in instances when stability 
considerations were crucial, and using a more accurate solver otherwise |21j . Sim- 
ilarly, we have since extended these methods to more general Monge- Ampere type 
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equations and have investigated the use of the various boundary conditions that 
arise in the mapping problem |22j . 
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